In [4]:
X.hist(figsize=(20, 15), bins=30)
plt.show()
from warnings import filterwarnings
from os import getcwd
from pathlib import Path
import numpy as np
import pandas as pd
import seaborn as sns
%matplotlib inline
import matplotlib.pyplot as plt
from sklearn.neighbors import LocalOutlierFactor
from factor_analysis import FactorAnalysis
from factor_analysis.plotting import scree_plot, create_loadings_heatmaps
from factor_analyzer.factor_analyzer import calculate_kmo
filterwarnings("ignore")
np.set_printoptions(suppress=True, precision=4)
pd.set_option("display.precision", 4)
Jede Beobachtung (Zeile) im Datensatz stellt eine sogenannte block group in Kalifornien dar.
Eine 'Block-Gruppe' ist eine statistische Aufteilung von des Volkszählungsamtes der USA, die zwischen 600 und 3000 Menschen umfassen sollten. Anstelle von Block-Gruppe werden wir hier meistens vereinfachend den Begriff Bezirk benutzen.
path = Path(getcwd(), "data", "cal_housing.data")
X = pd.read_csv(path, header=0, sep=",")
X.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 20640 entries, 0 to 20639 Data columns (total 9 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 longitude 20640 non-null float64 1 latitude 20640 non-null float64 2 housing_median_age 20640 non-null float64 3 total_rooms 20640 non-null float64 4 total_bedrooms 20640 non-null float64 5 population 20640 non-null float64 6 households 20640 non-null float64 7 median_income 20640 non-null float64 8 median_house_value 20640 non-null float64 dtypes: float64(9) memory usage: 1.4 MB
Insgesamt haben wir 20640 Beobachtungen (Bezirke) und 9 metrisch skalierte Merkmale.
Da wir für die Faktoranalyse auch metrisch skalierte Merkmale benötigen, haben wir in diesem Bereich keine Probleme.
X.describe()
| longitude | latitude | housing_median_age | total_rooms | total_bedrooms | population | households | median_income | median_house_value | |
|---|---|---|---|---|---|---|---|---|---|
| count | 20640.0000 | 20640.0000 | 20640.0000 | 20640.0000 | 20640.0000 | 20640.0000 | 20640.0000 | 20640.0000 | 20640.0000 |
| mean | -119.5697 | 35.6319 | 28.6395 | 2635.7631 | 537.8980 | 1425.4767 | 499.5397 | 3.8707 | 206855.8169 |
| std | 2.0035 | 2.1360 | 12.5856 | 2181.6153 | 421.2479 | 1132.4621 | 382.3298 | 1.8998 | 115395.6159 |
| min | -124.3500 | 32.5400 | 1.0000 | 2.0000 | 1.0000 | 3.0000 | 1.0000 | 0.4999 | 14999.0000 |
| 25% | -121.8000 | 33.9300 | 18.0000 | 1447.7500 | 295.0000 | 787.0000 | 280.0000 | 2.5634 | 119600.0000 |
| 50% | -118.4900 | 34.2600 | 29.0000 | 2127.0000 | 435.0000 | 1166.0000 | 409.0000 | 3.5348 | 179700.0000 |
| 75% | -118.0100 | 37.7100 | 37.0000 | 3148.0000 | 647.0000 | 1725.0000 | 605.0000 | 4.7432 | 264725.0000 |
| max | -114.3100 | 41.9500 | 52.0000 | 39320.0000 | 6445.0000 | 35682.0000 | 6082.0000 | 15.0001 | 500001.0000 |
Die Merkmale lassen sich also wie folgt beschreiben:
longitude: Längengrad eines Bezirks, Werte zwischen -124.35 und -114.31latitude: Breitengrad eines Bezirks, Werte zwischen -32.54 und 41.95housing_median_age: Median des Alters der Häuser in einem Bezirk, Werte zwischen 1 und 52 (Jahre)total_rooms: Gesamtzahl der Räume eines Bezirks, Werte zwischen 2 und 39320total_bedrooms: Gesamtzahl der Schlafzimmer eines Bezirks, Werte zwischen 1 und 6445population: Einwohnerzahl eines Bezirks, Werte zwischen 3 und 35682households: Gesamtzahl der Haushälte eines Bezirks, Werte zwischen 1 und 6082median_income: Median des Einkommens der Einwohner, Werte zwischen 0.5 und 15median_house_value: Median des Geldwertes der Häuser, Werte zwischen 15000 und 500 001 Dollar.Dabei fällt besonders der Wertebereich von median_income auf. Dies ist wahrscheinlich eine angepasste Skala und kein Einkommen in Dollar.
Wenn wir noch einmal die Beschreibung eines Bezirks anschauen, dann sehen wir, dass in diesem Datensatz
einige Ausreißer im Hinblick auf population vorliegen könnten. Denn ein Bezirk sollte hier zwischen 600 und 3000 Menschen umfassen.
Der Datensatz enthält keine fehlenden Werte, wie wir in der pd.info Methode sehen können.
Nun schauen wir uns einige univariate und bivariate Plots der metrischen Merkmale an, um ein Gefühl für die Daten zu bekommen.
X.hist(figsize=(20, 15), bins=30)
plt.show()
Hier sehen wir, dass die vier Merkmale total_rooms, total_bedrooms, population und households eine relativ hohe
Schiefe aufweisen.
Ebenso hat median_house_value bei circa 500 000 eine hohe Dichte. Dies könnte dadurch begründet sein,
dass der Wert 500 000 als obere Grenze benutzt wurde.
Bevor wir uns also bivariate Plots ansehen, werden wir mittels LocalOutlierFactor (LOF) versuchen, einige Ausreißer zu eliminieren.
labels = LocalOutlierFactor(n_neighbors=35).fit_predict(X)
outliers = X[labels == -1]
X = X[labels == 1]
print(f"Anzahl an mittels LOF als Ausreißer identifizierte Bezirke: {outliers.shape[0]}")
Anzahl an mittels LOF als Ausreißer identifizierte Bezirke: 701
X.describe()
| longitude | latitude | housing_median_age | total_rooms | total_bedrooms | population | households | median_income | median_house_value | |
|---|---|---|---|---|---|---|---|---|---|
| count | 19939.0000 | 19939.0000 | 19939.0000 | 19939.0000 | 19939.0000 | 19939.0000 | 19939.0000 | 19939.0000 | 19939.0000 |
| mean | -119.5843 | 35.6397 | 28.9650 | 2435.5189 | 499.6729 | 1325.3390 | 467.0708 | 3.8571 | 205807.2681 |
| std | 2.0012 | 2.1398 | 12.4808 | 1559.9463 | 318.4750 | 833.6129 | 292.7111 | 1.8869 | 114431.2155 |
| min | -124.3500 | 32.5400 | 1.0000 | 2.0000 | 1.0000 | 3.0000 | 1.0000 | 0.4999 | 36600.0000 |
| 25% | -121.8050 | 33.9300 | 19.0000 | 1435.0000 | 293.0000 | 781.0000 | 278.0000 | 2.5678 | 119750.0000 |
| 50% | -118.5000 | 34.2600 | 29.0000 | 2090.0000 | 427.0000 | 1149.0000 | 403.0000 | 3.5288 | 178600.0000 |
| 75% | -118.0200 | 37.7200 | 37.0000 | 3046.0000 | 624.0000 | 1668.0000 | 584.0000 | 4.7159 | 262500.0000 |
| max | -114.5500 | 41.9500 | 52.0000 | 16181.0000 | 3638.0000 | 8152.0000 | 3405.0000 | 15.0001 | 500001.0000 |
Wir können, sehen, dass der maximale Wert von population nun deutlich niedriger ist. Dennoch gibt es Bezirke mit sehr geringen Einwohnerzahlen (Minimum 3 Einwohner).
Schauen wir uns nun die Verteilungen der Merkmale mittels eines Pairplots genauer an.
sns.pairplot(X, diag_kind="kde", kind="scatter", plot_kws={"alpha": 0.3})
plt.show()
In diesem Pairplot sind univariate Plots auf der Hauptdiagonalen und bivariate Plots auf den nicht-diagonalen Elementen. Hier fallen besonders drei Gruppen auf:
Die Merkmale longitude und latitude weisen eine relativ hohe negative Korrelation auf.
Die Merkmale total_rooms, total_bedrooms, population und households korrelieren hoch miteinander. Dies ist jedoch
unter Beachtung der Bedeutung der Merkmale nicht überraschend, da eine hohe Anzahl an Menschen im Bezirk schließlich bedeuten muss, dass mehr
Räume existieren und insgesamt die Zahl der Haushälte wohl höher sein muss.
Das Merkmal housing_median_age könnte auch noch zu dieser Gruppe gezählt werden, da es eine leicht negative Korrelation zu diesen vier Merkmalen aufzeigt.
Die Merkmale median_income und median_house_value weisen ebenfalls eine sichtbare positive Korrelation auf, aber nicht so stark wie die in Gruppe 2.
Diese drei Gruppen sind unten noch einmal in gesonderten Pairplots dargestellt.
group1 = ["longitude", "latitude"]
group2 = ["total_rooms", "total_bedrooms", "population", "households"]
group3 = ["median_income", "median_house_value"]
sns.pairplot(X[group1], diag_kind="kde", kind="scatter", plot_kws={"alpha": 0.1}, size=4)
plt.show()
Hier sieht man, dass sich sehr viele Bezirke in zwei Regionen befinden. Dies erklärt die bimodale Struktur der beiden Merkmale longitude und latitude.
Weiter unten befindet sich ein Scatter-Plot, der die geographische Verteilung der Bezirke noch genauer zeigt.
sns.pairplot(X[group2], diag_kind="kde", kind="scatter", plot_kws={"alpha": 0.3}, size=3)
plt.show()
sns.pairplot(X[group3], diag_kind="kde", kind="scatter", plot_kws={"alpha": 0.3}, size=4)
plt.show()
Im letzten Pairplot fällt auch wieder die horizontale (bzw. vertikale) Linie bei median_house_value = 500 000 auf.
Zudem können wir uns noch die graphische Verteilung der Häuserpreise ansehen. Unten dargestellt ist ein Scatter-Plot, der die Lage der Bezirke in Abhängigkeit des Medians der Häuserpreise berücksichtigt.
fig, ax = plt.subplots(figsize=(15, 10))
sc = ax.scatter(x="longitude", y="latitude", c="median_house_value", cmap="viridis", data=X, alpha=0.4)
plt.colorbar(sc, label="median_house_value")
cities = ["Los Angeles", "San Francisco"]
xys = [(-118.243683, 34.052235), (-122.431297, 37.773972)]
xys_text = [(-121, 33.5), (-123, 35.5)]
for city, xy, xytext in zip(cities, xys, xys_text):
ax.annotate(city, xy=xy, xycoords='data',
xytext=xytext,
arrowprops=dict(facecolor='red', shrink=0.05),
horizontalalignment='right', verticalalignment='top',
)
plt.tight_layout()
plt.grid()
ax.set_xlabel("Längengrad")
ax.set_ylabel("Breitengrad")
plt.show()
Die Bezirke mit sehr teuren Häusern liegen am unteren Rand, also an der Küste. Zwei große Städte sind im Bild mit einem Pfeil markiert.
Zunächst müssen wir die Geeignetheit der Daten für die Faktoranalyse überprüfen. Dafür benutzen wir das Kaiser-Meyer-Olkin-Kriterium (KMO-Kriterium) für den kompletten Datensatz und das dazu verwandte Measure of sampling adequacy (MSA) für jedes einzelne Merkmal.
Generell wollen wir, dass der KMO-Wert über 0.5 ist und der MSA-Wert für jede Variable ebenfalls 0.5 nicht unterschreitet.
msa_values, kmo = calculate_kmo(X)
print(f"Der KMO-Wert beträgt {kmo:.4f}\n")
msa_df = pd.DataFrame(msa_values.reshape(-1, 1), index=X.columns, columns=["MSA"])
print(msa_df)
Der KMO-Wert beträgt 0.6537
MSA
longitude 0.4346
latitude 0.4325
housing_median_age 0.7372
total_rooms 0.8163
total_bedrooms 0.7182
population 0.8449
households 0.7390
median_income 0.3985
median_house_value 0.3872
Der KMO-Wert ist über 0.6, was auf eine akzeptable Qualität hinweist. Dies ist jedoch kein idealer Wert. Außerdem ist der MSA-Wert für 4 Variablen unter 0.5, jedoch ist der Wert für die anderen Merkmale gut.
Jetzt könnte man beispielweise die Merkmale mit einem MSA-Wert unter 0.5 entfernen. Für dieses Beispiel werden wir allerdings mit allen Variablen fortfahren.
Man sollte sich zudem noch die Korrelationsmatrix direkt anschauen. Dies tun wir im Folgenden mit einer Heatmap, da sie sich sehr gut als visuelle Repräsentation eignet.
fig = plt.figure(figsize=(10, 10))
corr = X.corr().round(2)
hm = sns.heatmap(data=corr, vmax=1, vmin=-1, cmap="RdBu_r", annot=True)
hm.set_xticklabels(X.columns, rotation=45)
hm.set_yticklabels(X.columns)
plt.tight_layout()
plt.show()
Wie wir schon im Pairplot gesehen haben, gibt es drei Gruppen von hoch korrelierten Merkmalen:
Die vier Merkmale in der Mitte, sowie housing_median_age mit einer leicht negativen Korrelation zu diesen vier Merkmalen.
Ebenso sind Longitude und Latitude negativ miteinander korreliert. Die Merkmale median_income und median_house_value sind positiv miteinander korreliert.
Dies deutet schon darauf hin, dass eine 3-Faktorlösung wahrscheinlich eine gute Wahl von k sein könnte. Im Folgenden werden dies genauer betrachten.
Um diese Frage zu beantworten, werden wir zunächst alle Faktoren mit der Hauptkomponentenmethode (Principal Component, PC) extrahieren und die Eigenwerte der Faktoren mithilfe eines 'Scree Plots' betrachten.
Dann können wir die Faktoren behalten, die einen Eigenwert größer Eins (Kaiser-Kriterium) haben, oder anhand eines 'Knicks' im Plot eine geeignete Faktorzahl erkennen.
fa = FactorAnalysis(n_factors=X.shape[1], method="pc").fit(X)
fig, ax = plt.subplots(figsize=(13, 8))
scree_plot(fa.eigenvalues_, ax)
ax.set_xlabel("Faktor")
ax.set_ylabel("Eigenwert")
plt.tight_layout()
plt.grid()
plt.show()
Hier können wir sehen, dass die ersten drei Faktoren einen Eigenwert größer als Eins haben. Dies macht im Hinblick auf die Korrelationsmatrix, die wir vorhin gesehen haben, wegen den drei 'Boxen' auch Sinn. Nach dem Kaiser-Kriterium sollten wir also ein 3-Faktor-Modell benutzen.
Der vierte Faktor ist jedoch nur minimal unter dem Eigenwert Eins, weshalb man diesen auch nicht direkt ausschließen sollte. Ein 'Knick' wäre bei $k = 5$ erkennbar. Wir werden also $k \in [3, 4, 5]$ ausprobieren und den 'Besten' auswählen.
Jetzt führen wir die eigentliche Faktorextraktion durch.
Dafür werden wir drei unterschiedliche Extraktionsmethoden miteinander vergleichen:
Dabei ist die letzte Variante wohl die am häufigsten eingesetzte Methode (unter diesen drei).
methods = [
("PC", FactorAnalysis(method="pc")),
("Nicht-iterierte PAF", FactorAnalysis(method="paf", max_iter=1)),
("Iterierte PAF", FactorAnalysis(method="paf", max_iter=50))
]
for n_factors in range(3, 6):
figsize = (10 + (1+n_factors)//2, 8)
create_loadings_heatmaps(X, methods, figsize, fa_params={"n_factors": n_factors})
plt.gcf().suptitle(f"Ladungsmatrizen eines {n_factors}-Faktor-Modells")
plt.show()
Hier sehen wir in jeder Zeile die Ladungsmatrizen der drei unterschiedlichen Methoden als Heatmap dargestellt. Wir stellen fest, dass die 3-Faktorlösung tatsächlich im Hinblick auf eine einfache Struktur eine gute Lösung darstellt. Die 4-Faktorlösung könnte jedoch auch noch eine valide Lösung sein. Lediglich die 5-Faktorlösung ist problematisch, da kein Merkmal hoch (absoluter Wert größer 0.4) auf den fünften Faktor lädt.
Ein Problem ist jedoch, dass das Merkmal housing_median_age nur bei der PC-Methode sehr hoch auf den vierten Faktor lädt
und bei einer 3-Faktorlösung nur eine moderate Ladung auf den ersten Faktor hat.
Dies deutet auf eine hohe spezifische Varianz hin, d.h. die Faktoren sind nicht gut in der Lage, die Varianz dieses Merkmals
zu erklären.
Wir werden also die 3-Faktorlösung genauer betrachten.
fa_params = {"n_factors": 3}
axes = create_loadings_heatmaps(X, methods, figsize=(10, 9), fa_params=fa_params, annotate=True)
plt.tight_layout()
plt.show()
Wir können sehen, dass bei der PC-Methode die Ladungen generell höher ausfallen und dass der zweite Faktor ein unterschiedliches Vorzeichen besitzt, im Gegensatz zu den anderen beiden Methoden.
Die iterierte und nicht-iterierte PAF-Methode sind sehr ähnlich zueinander. Jedoch ist die iterierte Variante oft in Hinblick auf die reproduzierte Korrelationsmatrix besser. Dies können wir anhand des root mean squared error (RMSE) untersuchen:
for method, fa in methods:
rmse = fa.get_rmse()
print(f"RMSE von {method}: {rmse:.4f}")
RMSE von PC: 0.0545 RMSE von Nicht-iterierte PAF: 0.0367 RMSE von Iterierte PAF: 0.0364
Der root mean squared error of residuals (RMSE) ist bei der iterierten PAF-Methode am geringsten, gefolgt von der nicht-iterativen Variante und der Hauptkomponentenmethode.
Schauen wir uns nun noch die Kommunalitäten der Variablen an. Dies zeigt uns zum Beispiel die print_summary Methode an:
# Zusammenfassung der iterierten PAF-Methode
methods[2][1].print_summary()
Fitted FactorAnalysis instance:
FactorAnalysis(heywood_handling='continue', initial_comm='smc',
is_corr_mtx=False, max_iter=50, method='paf', n_factors=3,
rotation=None)
Number of samples: 19939
Number of features: 9
Summary of estimated parameters
-------------------------------
F1 F2 F3 Communality Specific Variance
longitude 0.1382 -0.9356 -0.2010 0.9349 0.0651
latitude -0.1473 0.9598 0.0397 0.9445 0.0555
housing_median_age -0.3182 0.0033 0.0265 0.1020 0.8980
total_rooms 0.9430 0.0937 0.1561 0.9224 0.0776
total_bedrooms 0.9710 0.0946 -0.0662 0.9562 0.0438
population 0.8998 0.0242 -0.1189 0.8243 0.1757
households 0.9793 0.0903 -0.0444 0.9691 0.0309
median_income 0.0654 -0.1340 0.8297 0.7107 0.2893
median_house_value 0.0686 -0.1447 0.8055 0.6745 0.3255
Number of iterations: 7
Root mean squared error of residuals: 0.0364
Die spezifische Varianz von housing_median_age ist mit einem Wert von 0.8980 sehr hoch.
Das bedeutet, dass die Faktoren die Varianz dieses Merkmals gemeinsam nicht besonders gut erklären können.
Dies spiegelt sich ebenfalls in den geringen Ladungen auf die drei Faktoren wider.
Die spezifischen Varianzen bei den restlichen Merkmalen sind jedoch sehr niedrig,
was ein gutes Zeichen für die Qualität der Faktorlösung ist.
Da das Merkmal housing_median_age eine sehr hohe spezifische Varianz (geringe Kommunalität) aufweist
können wir auch ein 3-Faktor-Modell ohne diesem Merkmal anschauen.
Dieses kann den RMSE um knapp 44% reduzieren (im Vergleich sind die iterierten PAF-Methoden),
weshalb wir dieses Merkmal für die weitere Analyse entfernen.
X_without_age = X.drop(columns="housing_median_age", axis=1)
fa_without_age = FactorAnalysis(n_factors=3).fit(X_without_age)
perc = 1 - fa_without_age.get_rmse() / rmse
print(f"Der RMSE konnte durch Entfernen des Merkmals um {perc:.2%} reduziert werden")
Der RMSE konnte durch Entfernen des Merkmals um 43.50% reduziert werden
Den Unterschied zwischen den Methoden im RMSE können wir auch noch grafisch analysieren:
fig, axes = plt.subplots(ncols=2, figsize=(12, 6))
fitted_methods = [
("PC", FactorAnalysis(method="pc", n_factors=3).fit(X)),
("Iterierte PAF", FactorAnalysis(method="paf", n_factors=3, max_iter=50).fit(X))
]
for ax, (method, fa) in zip(axes, fitted_methods):
R = fa.corr_
R_hat = fa.get_reprod_corr()
abs_residuals = np.abs(R - R_hat)
mask = np.triu(np.ones_like(R))
ax.set_title(f"{method} (RMSE = {fa.get_rmse():.4f})", fontsize=11)
s = sns.heatmap(abs_residuals.round(2), cmap="BuGn", ax=ax, cbar=False, annot=True, square=True, mask=mask)
s.set_xticklabels(range(1, 10))
s.set_yticklabels(range(1, 10), rotation=0)
fig.suptitle("Residualmatrizen von zwei Extraktionsmethoden mit k=3")
fig.tight_layout()
plt.show()
Bevor wir zur Faktorrotation und -interpretation kommen, werden wir noch die verschiedenen initialen Schätzungen der Kommunalitäten in der (iterierten) PAF-Methode vergleichen. Interessant könnte dabei sein, ob die Wahl der initialen Schätzung einen Einfluss auf die finalen Ladungen hat.
paf_comparison_methods = [
("Nicht-iterierte PAF", FactorAnalysis(method="paf", max_iter=1)),
("PAF mit max_iter=3", FactorAnalysis(method="paf", max_iter=3)),
("Iterierte PAF", FactorAnalysis(method="paf", max_iter=50)),
]
figsize = (8, 6)
initial_communality_estimates = {
"smc": "Quadrierte multiple Korrelationen (SMC)",
"mac": "Maximale absolute Korrelationen (MAC)",
"ones": "Einsen"
}
for init_comm in initial_communality_estimates:
print(f"Initiale Schätzung: {initial_communality_estimates[init_comm]}")
create_loadings_heatmaps(X, paf_comparison_methods, figsize, fa_params={"n_factors": 3, "initial_comm" : init_comm})
plt.show()
print(f"Iterierte PAF hat {paf_comparison_methods[2][1].n_iter_} Iterationen benötigt.")
for method, fa in paf_comparison_methods:
print(f"RMSE von {method}: {fa.get_rmse():.4f}")
print("\n")
Initiale Schätzung: Quadrierte multiple Korrelationen (SMC)
Iterierte PAF hat 7 Iterationen benötigt. RMSE von Nicht-iterierte PAF: 0.0367 RMSE von PAF mit max_iter=3: 0.0364 RMSE von Iterierte PAF: 0.0364 Initiale Schätzung: Maximale absolute Korrelationen (MAC)
Iterierte PAF hat 3 Iterationen benötigt. RMSE von Nicht-iterierte PAF: 0.0373 RMSE von PAF mit max_iter=3: 0.0364 RMSE von Iterierte PAF: 0.0364 Initiale Schätzung: Einsen
Iterierte PAF hat 10 Iterationen benötigt. RMSE von Nicht-iterierte PAF: 0.0545 RMSE von PAF mit max_iter=3: 0.0370 RMSE von Iterierte PAF: 0.0364
Wir können sehen, dass nur geringe Unterschiede zwischen den unterschiedlichen initialen Schätzungen in den Ladungen feststellbar sind. Nur in der nicht-iterierten Variante der PAF-Methode können wir einige Unterschiede, vor allem im zweiten Faktor feststellen. Beispielsweise hat der zweite Faktor hier ein unterschiedliches Vorzeichen, jedoch nur, wenn Einsen als Kommunalitätsschätzung benutzt wurden. In diesem Fall ist das Ergebnis identisch zur Hauptkomponentenmethode.
Wir stellen fest, dass die iterierte Variante jedoch eine unterschiedliche Anzahl an Iterationen benötigt, bis das Konvergenzkriterium erreicht wird. Am langsamsten ist es auf diesem Datensatz im Falle von Einsen (10 Iterationen) und am schnellsten ist es bei den maximalen absoluten Korrelationen (MAC) (nur 3 Iterationen). Ist das Konvergenzkriterium erfüllt, sind die Ladungen jedoch bei allen drei initialen Schätzungen weitestgehend identisch.
Jetzt rotieren wir die Ladungen mit der Varimax-Methode und versuchen, die Faktoren zu interpretieren.
methods = [
("Unrotiert", FactorAnalysis(method="paf", rotation=None)),
("Varimax-Rotation", FactorAnalysis(method="paf", rotation="varimax"))
]
fa_params = {"n_factors": 3}
fig = create_loadings_heatmaps(X_without_age, methods, fa_params=fa_params)
plt.tight_layout()
Müssten wir die Faktoren interpretieren, könnte man sagen, dass
Als letzten Schritt können wir noch einen Blick auf die geschätzten Faktorwerte werfen. Wie hätte als der Bezirk $i$ die drei oben beschriebenen Faktoren bewertet?
Die hier benutzte Methode zur Schätzung der Faktorwerte ist die Regressionsmethode, welche eine multivariate lineare Regression benutzt, um die Faktorwerte zu schätzen.
scores = FactorAnalysis(n_factors=3).fit_transform(X_without_age)
scores = pd.DataFrame(scores, columns=["Größe", "Standort", "Wohlstand"])
scores.head()
| Größe | Standort | Wohlstand | |
|---|---|---|---|
| 0 | -1.1775 | 0.7685 | 2.2611 |
| 1 | -0.9650 | 0.8176 | 1.7476 |
| 2 | -0.9210 | 0.9032 | 1.2176 |
| 3 | -0.8342 | 0.9682 | 0.8220 |
| 4 | -1.0621 | 0.9492 | 0.5564 |
scores.std(axis=0)
Größe 0.9954 Standort 0.9857 Wohlstand 0.9195 dtype: float64
Hier weisen die Faktorwerte keine Einheitsvarianz beziehungsweise Standardabweichung von eins auf, weil die PAF-Methode mit quadrierten multiplen Quadraten als initiale Schätzung verwendet wurde.
Benutzt man jedoch die Hauptkomponentenmethode, so weisen die Faktorwerte eine Standardabweichung von Eins auf.
scores = FactorAnalysis(method="pc", n_factors=3).fit_transform(X_without_age)
scores.std(axis=0)
array([1., 1., 1.])